%print the result of question2
clc; clear ; close all;
a = 0;
b = 1;
T = 2;
cfl = 0.2;
u01 = @(x) sin(2*pi*x);
u02 = @(x) (x >= 1/4 & x <= 3/4);
h = [1/8, 1/16, 1/32, 1/64, 1/128];
k = [1,2];
errors_inf_RK2 = zeros(2,5);
errors_inf_RK3 = zeros(2,5);
errors_L2_RK2 = zeros(2,5);
errors_L2_RK3 = zeros(2,5);
[b_fun, db_fun] = pro_basis();
%创建一个矩阵用来存h=1/128时，间断初值的数值解
tmp_x = a : h(5) : b;
u_discon_num = zeros(4,length(tmp_x)-1);
x_half = (tmp_x(1:end-1) + tmp_x(2:end))/2;
discon_id = 1; 
u_discon_true = u02(x_half);
%solve the first question
for i = 1 : 2
    kk = k(i);
    for q = 1 : 5
        A01_coef_RK2 = DG_CL_generalize(u01, kk, h(q), T, cfl,a, b, 0);
        A01_coef_RK3 = DG_CL_generalize(u01, kk, h(q), T, cfl, a, b, 1);
        X = a : h(q) : b;
        x_half = (X(1:end-1) + X(2:end))/2;
        N = length(X) - 1;
        if q == 5
            A02_coef_RK2 = DG_CL_generalize(u02, kk, h(q), T, cfl,a, b, 0);
            A02_coef_RK3 = DG_CL_generalize(u02, kk, h(q), T, cfl, a, b, 1);
            for cell = 1 : N
                x_mid = 0.5*(X(i) + X(i+1));
                for l = 1 : kk+1
                    u_discon_num(discon_id,cell) = u_discon_num(discon_id,cell) + A02_coef_RK2(l,cell)*b_fun{l}(0);
                    u_discon_num(discon_id+1, cell) = u_discon_num(discon_id+1, cell) + A02_coef_RK3(l,cell)*b_fun{l}(0);
                end
            end
            discon_id = discon_id + 2;
        end
        u1 = @(x) u01((x-T));
        for cell = 1 : N
            x_left = X(cell);
            x_right = X(cell+1);
            x_mid = (x_left + x_right)/2;
            u_num_RK3 = @(x) 0.*x;
            u_num_RK2 = @(x) 0.*x;
            for l = 1 : kk+1
                Pl = @(x) b_fun{l}((2*(x - x_mid)/h(q)));
                u_num_RK3 = @(x) u_num_RK3(x) + A01_coef_RK3(l,cell)*Pl(x);
                u_num_RK2 = @(x) u_num_RK2(x) + A01_coef_RK2(l,cell)*Pl(x);
            end
            error_RK3 = @(x) (-1)*abs(u_num_RK3(x) - u1(x));
            error_RK2 = @(x) (-1)*abs(u_num_RK2(x) - u1(x));
            tmp_RK3 = fminbnd(error_RK3, x_left, x_right);
            tmp_RK3 = -error_RK3(tmp_RK3);
            if tmp_RK3 > errors_inf_RK3(i,q)
                errors_inf_RK3(i,q) = tmp_RK3;
            end
            tmp_RK2 = fminbnd(error_RK2, x_left, x_right);
            tmp_RK2 = -error_RK2(tmp_RK2);
            if tmp_RK2 > errors_inf_RK2(i,q)
                errors_inf_RK2(i,q) = tmp_RK2;
            end
            errors_L2_RK2(i,q) = errors_L2_RK2(i,q) + quadgk(@(x) error_RK2(x).^2, x_left,x_right);
            errors_L2_RK3(i,q) = errors_L2_RK3(i,q) + quadgk(@(x) error_RK3(x).^2, x_left,x_right);
        end
        errors_L2_RK2(i,q) = sqrt(errors_L2_RK2(i,q));
        errors_L2_RK3(i,q) = sqrt(errors_L2_RK3(i,q));
    end
end
orders_inf_RK3 = log2(errors_inf_RK3(:,1:end-1)./errors_inf_RK3(:,2:end));
orders_inf_RK2 = log2(errors_inf_RK2(:,1:end-1)./errors_inf_RK2(:,2:end));
orders_L2_RK3 = log2(errors_L2_RK3(:,1:end-1)./errors_L2_RK3(:,2:end));
orders_L2_RK2 = log2(errors_L2_RK2(1,1:end-1)./errors_L2_RK2(1,2:end));
%输出误差和阶数
fprintf("k = 1, RK2\n");
fprintf("L2 误差：");
for i = 1:length(errors_L2_RK2(1,:))
    fprintf(" %e", errors_L2_RK2(1,i));
end
fprintf("\n");

fprintf("L2 阶数：%s\n", num2str(orders_L2_RK2));
fprintf("inf 误差：");
for i = 1:length(errors_inf_RK2(1,:))
    fprintf(" %e", errors_inf_RK2(1,i));
end
fprintf("\n");

fprintf("inf 阶数：%s\n", num2str(orders_inf_RK2(1,:)));

fprintf("\nk = 1, RK3\n");
fprintf("L2 误差：");
for i = 1:length(errors_L2_RK3(1,:))
    fprintf(" %e", errors_L2_RK3(1,i));
end
fprintf("\n");

fprintf("L2 阶数：%s\n", num2str(orders_L2_RK3(1,:)));
fprintf("inf 误差：");
for i = 1:length(errors_inf_RK3(1,:))
    fprintf(" %e", errors_inf_RK3(1,i));
end
fprintf("\n");

fprintf("inf 阶数：%s\n", num2str(orders_inf_RK3(1,:)));

fprintf("\nk = 2, RK2\n");
fprintf("L2 误差：");
for i = 1:length(errors_L2_RK2(2,:))
    fprintf(" %e", errors_L2_RK2(2,i));
end
fprintf("\n");
fprintf("inf 误差：");
for i = 1:length(errors_inf_RK2(2,:))
    fprintf(" %e", errors_inf_RK2(2,i));
end
fprintf("\n");

fprintf("\nk = 2, RK3\n");
fprintf("L2 误差：");
for i = 1:length(errors_L2_RK3(2,:))
    fprintf(" %e", errors_L2_RK3(2,i));
end
fprintf("\n");

fprintf("L2 阶数：%s\n", num2str(orders_L2_RK3(2,:)));
fprintf("inf 误差：");
for i = 1:length(errors_inf_RK3(2,:))
    fprintf(" %e", errors_inf_RK3(2,i));
end
fprintf("\n");

fprintf("inf 阶数：%s\n", num2str(orders_inf_RK3(2,:)));
%输出间断初值情形的图像
subplot(2,2,1);
plot(x_half,u_discon_true,"-b",x_half,u_discon_num(1,:),":r");
title("k = 1, RK2");
subplot(2,2,2);
plot(x_half,u_discon_true,"-b",x_half,u_discon_num(2,:),":r");
title("k = 1, RK3");
subplot(2,2,3);
plot(x_half,u_discon_true,"-b",x_half,u_discon_num(3,:),":r");
title("k = 2, RK2");
subplot(2,2,4);
plot(x_half,u_discon_true,"-b",x_half,u_discon_num(4,:),":r");
title("k = 2, RK3");
%plot(x_half,u_true, x_half, u_num1);